This example uses the rayshader package to visualize the landscape and landform surrounding a pedon, while also displaying a 3d representation of the pedon horizonation, depths and colors.

Using the latitude and longitude we create a buffer of 2 kilometers to capture a bounding box of the surrounding area. This is then used as an input for querying the elevation data from an online web service. The extent is then cropped to limit the area of interest. After the elevation data is obtained, the shading can be generated.

First create a connection to the local NASIS database. You could also use fetchNASIS() to obtain the site data.

library(DBI)
con <- dbConnect(odbc::odbc(), "nasis_local", timeout = 10, uid = "NASISSQLRO", 
    pwd = "nasisRe@d0n1y365")
sers <- params$series

Although any site/pedon could be used, in this case, the Official Series Description for the Marshall series and its associated site will be used. A custom SQL query was created to obtain the latitude and longitude of the OSD pedon type with a taxon name = Marshall from NASIS.


SELECT TOP 1 site.latstddecimaldegrees, site.longstddecimaldegrees, site.siteiid --if more than one exists, this insures only one is selected
FROM site
JOIN siteobs ON site.siteiid = siteobs.siteiidref
JOIN pedon ON pedon.siteobsiidref = siteobs.siteobsiid
WHERE pedontype = 7 -- osd pedon type
AND taxonname = ?sers

Required libraries:

library(soilDB)
library(aqp)
library(sf)
library(raster)
library(dplyr)
library(rayshader)
library(rgl)
library(elevatr)

Two custom functions were developed to convert and plot a soil profile collection object into a 3D representation in rgl. The first function creates the wireframe that helps to distinguish horizon boundaries. The second function builds the faces of the pedon object and shades it with the appropriate colors.

source("build3d_profwhex.R")
source("build3d_profhex.R")

#make an sf object from the lat/long
s <- st_as_sf(nq, coords =c( "longstddecimaldegrees", "latstddecimaldegrees"), crs = 4326)

#create a 2 km buffer around the point
se <- st_buffer(s, 2000)

#use the buffer area to query the elevation data
r11r <- get_elev_raster(se, z = 14) # different zoom levels here will change the results and formatting significantly

#crop raster to match the extent of the buffer zone, (which centers the raster at the lat/long of the pedon) then crop again to set specific dimensions
cc <- crop(r11r, st_bbox(se))
cc <- crop(cc, extent(cc, ((nrow(cc)/2)-250),((nrow(cc)/2)+250), ((ncol(cc)/2)-250), ((ncol(cc)/2)+250) ))

#convert the raster to a matrix so it can be used in rayshader
ccmat <- raster_to_matrix(cc)

#compute ambient shade for the matrix
ambmat <- ambient_shade(ccmat, sunbreaks = 24, maxsearch = 300, multicore=TRUE)

n_frames <- 180

# thetavalues <- -90 + 45 * cos(seq(0, 2*pi, length.out = n_frames))
thetavalues <- -90 + 180 * cos(seq(0, pi, length.out = n_frames))

img_frames <- paste0(params$series, "_series", seq_len(n_frames), ".png")
for (i in seq_len(n_frames)) {
  message(paste(" - image", i, "of", n_frames))
  rgl::clear3d()  
  w <- build3dprofilewhex(params$series, cc)
  
render_points(extent = extent(cc), heightmap = ccmat, lat = nq$latstddecimaldegrees, long = nq$longstddecimaldegrees, zscale = 1, color = "red")

tst <- ccmat %>%
  sphere_shade(zscale = 10, texture = "imhof1") %>%
  add_shadow(ray_shade(ccmat, zscale = 10, maxsearch = 300)) %>%
  add_shadow(ambmat, max_darken = .8) %>%
  plot_3d(ccmat, zscale = 1, fov = 0, theta = thetavalues[i], zoom = .9, phi = 45, windowsize = c(768, 768), baseshape = "hex")
# render_label(ccmat, x = label$pos$x, y = label$pos$y, z = (cc@data@max/2.5)+200, zscale = 2.5, text = label$text, textsize =1, linewidth = 2, freetype = FALSE )

# axis3d(edge = "y++", color = "black", pos = c(90,NA,-60), labels = fh1$top, at = (max(fh1$bottom)-fh1$top)/2.5 + cc@data@max/2.5 +(250 - max(fh1$bottom))/2.5)

# axis3d(edge = "y++", color = "black", pos = c(90,NA,-60), labels = fh1$bottom, at = (max(fh1$bottom)-fh1$bottom)/2.5 + cc@data@max/2.5 +(250 - max(fh1$bottom))/2.5)

# axis3d(edge = "y--", color = "black", pos = c(60,NA,-90), labels = fh1$hzname, at = ((max(fh1$bottom)-((fh1$bottom+fh1$top)/2)) /2.5) + cc@data@max/2.5 +(250 - max(fh1$bottom))/2.5)



t <- build3dprofilehex(params$series, cc)

  render_snapshot(img_frames[i])
  rgl::clear3d()
}

magick::image_write_gif(magick::image_read(img_frames), 
                        path = paste0(params$series, "_series_ani_hex.gif"), 
                        delay = 12/n_frames)
## [1] "M:\\GitHub\\rayshader-examples\\OSD-pedon-loc-viz\\Marshall_series_ani_hex.gif"